Itinerant ferromagnetism in a Fermi gas with contact interaction: 
Magnetic properties in a dilute Hubbard model 
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Ground state properties of the repulsive Hubbard model on a cubic lattice are investigated by 
means of the auxiliary-field quantum Monte Carlo method. We focus on low-density systems with 
varying on-site interaction U/t, as a model relevant to recent experiments on itinerant ferromag- 
netism in a dilute Fermi gas with contact interaction. Twist-average boundary conditions are used 
to eliminate open-shell effects and large lattice sizes are studied to reduce finite-size effects. The 
sign problem is controlled by a generalized constrained path approximation. We find no ferromag- 
netic phase transition in this model. The ground-state correlations are consistent with those of a 
paramagnetic Fermi liquid. 
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Introduction - The study of ferromagnetism has a long 
history in physics. At the microscopic level, the forma- 
tion of ferromagnetic order is a consequence of strong in- 
teractions. Heisenberg first pointed out that an exchange 
interaction that lowers the energy of a pair of parallel 
spins would favor ferromagnetism. However a localized- 
spin mechanism cannot be fully responsible for ferromag- 
netism in transition metals, for instance iron and nickel, 
where electrons are extended. Both the interactions and 
the delocalized nature of electrons have to be taken into 
account at a more fundamental level. 

A generic description of itinerant ferromagnetism is 
given by the three-dimensional (3-D) Hubbard model [![: 
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The operator c\ a (ci a ) creates (annihilates) an electron 
with spin a (a =j, 1 ), i enumerates the sites in an 
N = L 3 lattice, and (ij) denotes a sum of nearest neigh- 
bor pairs. The parameter t is the nearest-neighbor hop- 
ping amplitude, and U > is the on-site interaction 
strength. The total density is n = (JVf + N±)/N. This 
simple Hamiltonian contains both the itinerant character 
and local repulsion. However, because neither of the two 
terms alone favors ferromagnetic ordering, the magnetic 
correlations in the Hubbard model is not obvious. 

The first evidence of ferromagnetism in the Hubbard 
model was discussed by Nagaoka [|[ and by Thouless Q . 
They showed that the ground state with a single hole 
in any finite bipartite lattice with U — > oo (and periodic 
boundary conditions) is fully polarized. Subsequent stud- 
ies indicate that the stability of the state with more holes 
can depend on system size [4j, and boundary conditions 
Q . The critical dopingfor the onset of ferromagnetism is 
still an open question [J, [|| 0] • Away from infinite- U, the 
existence of ferromagnetism at non-zero density is less 
certain. Whether ferromagnetism is a generic property 
of the Hubbard model is still not answered. 



Rapid experimental progress in cold atoms has opened 
a new avenue for exploring the physics of itinerant ferro- 
magnetism. In a recent experiment aimed to simulate the 
Stoner Hamiltonian (i.e. spin 1/2 fermions in continuous 
space interacting with a repulsive contact potential), a 
dilute gas of two hyperfine states of 6 Li atoms are tuned 
to interact via large positive scattering lengths. Signa- 
tures of ferromagnetic instability [|| have generated a lot 
of theoretical interest l9T-[l2ll. 



The Hubbard Hamiltonian in Eq. (fTJ) gives a reasonable 
representation of the Stoner Hamiltonian on a lattice. As 
the density n — > it seems clear that no ferromagnetism 
exists in the model in Eq. (fT]) [13], since the maximum 
scattering length is bounded by ~ 1/3.173 lattice spac- 
ing [lil . Il5j . However, at low but not zero density, the 
magnetic properties and the phase diagram of the 3-D 
Hubbard model are not clear. We address this question 
here by very accurate many-body simulations with the 
constrained path Monte Carlo (CPMC) method. 

Method - The CPMC method [ll-QJl projects the 
many-body ground state \^o) from a trial wave func- 
tion \^>t) by repeated application of an imaginary-time 
propagator e~ ArH (At is the Trotter time step), pro- 
vided that \^t) satisfies (^oI^t) 0. The propagator is 
decomposed into e- ArH « e -ArH 1 /2 e -ArH 2e -ArH 1 /2 + 

C(Ar 3 ), where Hi and Hi are one- and two-body parts 
of H respectively. The two-body part e - ArH 2 is fur- 
ther decou pled into a sum over one-body projectors in 
Ising fields 19]. This leads to a formally exact expression 
e~ ArH = J2{ x } -P({ x })^({ x })> where {x} is a collection 
of N Ising fields, P({x}) is their probability distribution, 
and -B({x}) is a one-body projector. The multidimen- 
sional summation is carried out efficiently by importance- 
sampled random walks with non-orthogonal Slater deter- 
minants (SDs), where the one-body projectors £?({x}) 
propagate one SD into another. 

The fermion sign problem is controlled approximately 
by the constrained path approximation. (16j The many- 
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body ground state is given by |^o) = S</> ^Wl^): where 
|0) are SDs sampled by the QMC, and their probabil- 
ity distribution determines the weight factors w(<f>). Be- 
cause the Schrddinger equation is linear, l^o) is degen- 
erate with — l^o)- In a random walk, the SDs can move 
back and forth between the two sets of solutions. The ap- 
pearance of the two sets with opposite signs in the Monte 
Carlo samples is the origin of the sign problem. To con- 
trol the problem, the walker is required to satisfy the 
constraint (^tI^) > in the course of the random walk. 
This is the only approximation in our method. More for- 
mal discussions of the theoretical basis of the generalized 
constrained path approximation and benchmarks can be 
found elsewhere [ljl [l7| . In the Hubbard model, the en- 
ergy at U = At is typically within < 0.5% of the exact 
diagonalization result [1 81 ] . Extensive benchmarks of this 
approach for molecules and solids are in Refs. |20l l2lj. 

The constrained path approximation is similar in spirit 
to the fixed-node approximation in the diffusion Monte 
Carlo (DMC) method [H[23j], which has been used for 
all recent simulation work on the problem of itinerant 
ferromagnetism in the Stoner model pil HH, H3| ■ In fixed- 
node DMC one uses a real-space trial function ^>p(R) to 
determine the sign of the ground-state wave function. 
The random walks, which involve movements of electron 
coordinates R (a 3 (N^ + iVjJ-dimensional vector), are 
constrained to the region where <3/t(-R) > 0. Since, in 
CPMC the random walks take place in the space of SDs, 
where fermionic statistics are automatically maintained, 
the sign problem is reduced. As a result, the constrained 
path approximation is less sensitive to | V^t) and typically 
has smaller systematic errors. 

In this work we apply twist-averaged boundary con- 
ditions (TABCs) [25(. Under TABCs, the wave func- 
tion gains a phase when electrons wind around the 
periodic boundary conditions: ^(...,rj + L, . . .) = 

e* L $(. . . . . .), where L is the unit vector along L, 

and © = {0 x ,9 y ,9 z ) are random twists over which we 
average. A simple generalization of the CPMC method 
can be made to handle the overall phase that arises from 
TABC [HI EH- As an additional benchmark for the 
present work, we studied several low-density L = A sys- 
tems in detail. For example, at U = 16t with n = 0.25, 
the CPMC energy, averaged over 1000 ©-points, agrees 
to better than 0.2% with exact diagonalization. 



Energy - We first compare the ground-state energy of 
an unpolarized system (iV-j- = N±) with that of a fully 
polarized state at the same total density. The results 
are summarized in Fig. [1] Because electrons of the same 
spin do not interact, the energy of the fully polarized 
state, 6fmi is purely kinetic and does not depend on 
U . In mean- field (MF) theory, the energy of a system 
with n a — N a /N (with n = rif + nj.) is eMF(U,n) = 
[eo(iVf)iVf + e (ni)ni + £/n t nJ/n, where e (n a ) is the 
energy of the fully polarized system at density n a . At 
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FIG. 1. (color online). Ground state energy per particle e as 
a function of interaction strength U/t at n — 0.25 (left) and 
n = 0.0625 (right). Symbols represent ecPMC- Dashed (blue) 
line corresponds to the energy of a saturated ferromagnetic 
state (efif). eMF energy is represented by the thick solid 
(green) line, ep (perturbation theory [26j]) is plotted by dot- 
dashed line. 



n = 0.25, MF predicts a paramagnetic to ferromagnetic 
phase transition at U = 13. 9^. This is to be compared 
to the corresponding transition point kpa ~ 7r/2 in the 
continuum Stoner Hamiltonian, where kp — {iir 2 n) 1 ^ 
is the Fermi wave vector, and a is the scattering length 
in continuum. When the system density is lowered to 
n = 0.0625, the MF transition in the Hubbard model 
is at a larger interaction, U = 29. 3t. The equation of 
state has also been obtained from perturbation theory 
for an unpolarized system [26[: ep(U,n) — eMF(U,n) + 
e c (U, n). The last term, e c (U, n), is the correlation energy 
estimated to 0(U 2 ). The result is also included in Fig.[T] 
The CPMC result for the ground state energy ec pmc 
is obtained by averaging over twist-angles. The energies 
calculated from different lattice sizes are shown by differ- 
ent symbols in Fig. [T] It can be seen that our remaining 
finite-size errors are negligible on this scale. Free-electron 
trial wave functions are used for the constraint. In a few 
cases we have also checked with unrestricted Hartrcc- 
Fock trial wave functions, which gave statistically indis- 
tinguishable CPMC energies. The energies shown are for 
finite time steps, with At satisfying UAt < 0.2. The 
residual Trotter error is O(10~ 2 ), smaller than the sym- 
bol size. 

We see that MF theory, which gives a reasonable es- 
timate of the energy at small U, quickly shows severe 
deviations as the interaction becomes stronger. The per- 
turbation result, ep{U, n), gives an improved estimate of 
energy for small U, but deviates once the system enters 
the intermediate interaction regime U > 5t. At the MF 
transition point, the CPMC energy is significantly lower 
than epM- Indeed the CPMC energy remains lower than 
epM across the entire range of U simulated. No indica- 
tion of a ferromagnetic transition is seen. 

Individual components of the energy are shown in 
Fig. [21 As U increases, electrons in the unpolarized sys- 
tem occupy higher momentum states, outside the Fermi 
level, which increases the kinetic energy compared to the 
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FIG. 2. (color online). Kinetic (left panel) and interaction 
(right panel) energies as a function of interaction strength at 
density n — 0.0625. Symbols are the CPMC data obtained 
on an 8 3 lattice. Lines are defined in the same way as in 
Fig. [1] The inset on the right shows the double occupancy, 
normalized to 1 at U = 0. 



MF result. This enables the system to drastically de- 
crease the interaction energy, by suppressing double oc- 
cupancy. The net effect is that the total energy is greatly 
reduced and remains below euF and epu- 



Correlation junction - To probe the nature of the 
ground state, we examine the spin-dependent pair cor- 
relation function: 



flW' ( r ) 
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The CPMC expectations are evaluated by the back- 
propagation technique [l|| [27j ■ We average over different 
r's to obtain g a a'(r), with r = |r| since the correlation 
function is primarily a function of distance in the para- 
magnetic or ferromagnetic phases. 



The anti-parallel pair correlation g^(r) is a constant 
in a non-interacting system or in the MF solution. In 
the presence of interaction, a correlation hole is created 
surrounding each electron. At n = 0.0625, the size of 
the correlation hole is r cor < y/3. As U is increased, 
the correlation hole becomes deeper, as illustrated in the 
left panel in Fig. [3] Compared to gf±(r), the change in 
the parallel-spin pair correlation g^f(r) is less dramatic 
from the MF or non- interacting result. Strong interaction 
does appear to increase gff (r) slightly at short distance. 
However, the correlation remains much less than that in 
the FM case. 

Momentum distribution - The creation of correlation 
hole is a result of minimizing the interaction energy. Elec- 
trons of opposite spins rearrange their relative positions 
to reduce the potential energy. The cost of the rear- 
rangement is the kinetic energy increase, as discussed 
earlier. This can also be observed in Fig. [U where the 
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FIG. 3. (color online). Left: Anti-parallel spin-spin pair 
correlation function of an unpolarized state at different in- 
teraction strengths. Right: Comparison of the parallel pair 
correlation functions of the fully polarized state (FM) and 
the unpolarized states at different interaction strengths. The 

inset shows Aptf(r) = <7-|-1~( r ) ~ 5?t( r )' wnere 5?t( r ) ls the 
correlation function of the unpolarized non-interacting sys- 
tem. In both panels, the system is an 8 3 lattice at density 
n = 0.0625. 



momentum distribution is shown for different inter- 
action strengths. We have plotted rik as a function of 
the single-particle energy level e(k) = 2^ Q=2 , y2 [l — 
cos(fc Q + Q a /L)j, in units of the Fermi energy ep. Each 
curve contains the result of from multiple 0-points. 
At U = 4t, the distribution is very close to the non- 
interacting momentum distribution with only a few low 
lying excitations near the Fermi surface (FS). As U is in- 
creased, more higher k states are populated outside the 
FS. In 7ik a jump appears at ep which can be read off 
directly in our finite size simulations. The jump indicates 
that the system is a normal Fermi liquid, with the value 
of the jump proportional to the renormalization factor Z . 
Its precise value can be determined with more extensive 
calculations and finite size scaling. 

Discussion - Although we have focused on the ground 
state of a homogeneous Fermi gas, it is not difficult to 
extend the results to the case with an external trap. For 
example, the kinetic energy results in Fig. [^indicate that, 
with a trap, there would be a minimum in the curve of 
ex versus interaction strength, as observed in the ex- 
periment [H (see also discussion in Ref. 9]). Effects of 
confinement on the kinetic energy have been investigated 
in detail by CPMC for trapped Bose gases [Hj]. The MF 
kinetic energy was shown to decrease monotonically be- 
cause the gas expands in the trap as the scattering length 
a is increased; on the other hand, correlation effects lead 
to an increase of e^, similar to Fig. [2j This competition 
results in a non-monotonic curve, with a minimum in the 
kinetic energy at a finite scattering length. 



The Hubbard model, of interest in its own right, con- 
tains some of the same features (namely itinerant elec- 
trons and local interaction) as the continuum Stoner 
Hamiltonian. However, there are differences with respect 
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FIG. 4. The momentum distribution rik for values of U plot- 
ted as a function of the single particle energy. Two lattice 
sizes are shown, at the density n — 0.0625. In each system, 
we average over 10 random twist angles. 

to the experiment worth emphasizing. The experiment 
is in the continuum, using an attractive (negative) inter- 
action with an effective positive scattering length for the 
excited state that describes the prepared state. (How- 
ever, there are questions whether such an effective de- 
scription is appropriate 12].) In our simulation, we use 
a discretized representation, with positive on-site inter- 
action. As mentioned before, the lattice model leads to 
a scattering length bounded by roughly the lattice spac- 
ing. Using the above values for the maximum scattering 
length and kp in the unpolarized phase, we find that 
kpa < 1.0371 1 / 3 . (For the transition to the ferromag- 
netic phase, n < 1 since one cannot have two like-spin 
fermions on a single site. To focus on the dilute limit, we 
have done calculations for up to n = 0.5.) 

Recently the problem of itinerant ferromagnetism in 
repulsive Fermi gases has been studied by several groups 
[IS ED, HI using the DMC method with the fixed-node 
approximation. These calculations all found the exis- 
tence of a ferromagnetic instability. The DMC calcula- 
tions were all done in the continuum, while the present 
calculation is for the Hubbard (lattice) model. In the 
DMC simulations, the atomic interaction is modeled by a 
repulsive potential whose range is determined by the scat- 
tering length. We note that since the scattering length 



diverges near resonance, the range of potential (or the 
range of the node in the Jastrow when a negative interac- 
tion is used) can become very large. Note that the hard- 
sphere-likc interaction is only between unlike spins. As 
the scattering length approaches the interparticle spac- 
ing, there is a strong tendency to separate into a spin-up 
and spin-down domains, to lower the interaction energy; 
i.e. it favors ferromagnetism. 

Of lesser importance, the DMC fixed-node errors in the 
calculated ground state energies bias the result in favor 
of ferromagnetism, since nodal surfaces for the ferromag- 
netic state are more accurate than the spin unpolarized 
state [28[ . Although the constrained path error from our 
calculations could also be biased, previous calculations 
indicate 

[3 [U [Hi that the systematic error in CPMC 
is smaller than the fixed-node error from single determi- 
nant trial wave functions used in these calculations. 

Summary - We have examined the magnetic proper- 
ties in the ground state of the dilute 3D Hubbard model, 
using the CPMC method and twist-averaged boundary 
conditions. Our simulation results indicate that there 
is no ferromagnetic instability in this model with strong 
on-site repulsions for densities up to 0.5. The ground 
state appears to be a paramagnetic Fermi liquid. The 
total energy is effectively lowered by electron correlation 
which, while increasing the kinetic energy, can strongly 
suppress double occupancy to lower the interaction en- 
ergy. In the presence of a trap, the kinetic energy can be 
decreased by the expansion of the gas due to repulsive 
interaction. A kinetic energy minimum, which was ob- 
served in the experiment, can be understood in terms of 
the competitions between these effects. We have also dis- 
cussed the difference between our calculations and recent 
results from DMC simulations, as well as connections and 
differences with the Fermi gas itinerant ferromagnetism 
experiments. 
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